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ABSTRACT 

We compute the production rate of photons in the ionizing Lyman continua (LyC) of H I (A < 
912 A), He I (A < 504 A), and He II (A < 228 A) using recent stellar evolutionary tracks coupled to a 
grid of non-LTE, line-blanketed WM-basic model atmospheres. The median LyC production efficiency 
is QcyC = (6± 2) X 10®° LyC photons per Mq of star formation (range [3.1 — 9.4] x 10®°) corresponding 
to a revised calibration of 10®° °=*=°'^ photons s“^ per Mq yr“^. Efficiencies in the helium continua 
are Qnei ~ 10®° photons Mq^ and Queii ~ 10®® photons at solar metallicity and larger at low 
metallicity. The critical star formation rate needed to maintain reionization against recombinations 
at z ss 7 is psFR = (0.012 Mq yr“^ Mpc“°)[(l -I- z)/8]°(C'H/3)(0.2//esc) for fiducial values of IGM 
clumping factor Ch ~ 3 and LyC escape fraction /esc ~ 0.2. The boost in LyC production efficiency is 
an important ingredient, together with metallicity, Ch, and /esc, in assessing whether IGM reionization 
was complete by z « 7. Monte-Garlo sampled spectra of coeval starbursts during the first 5 Myr have 
intrinsic flux ratios of E'(1500)/F(900) ss 0.4 — 0.5 and F(912“)/F(912+) « 0.4 — 0.7 in the far-UV 
(1500 A), the LyC (900 A), and at the Lyman edge (912 A). These ratios can be used to calibrate the 
LyC escape fractions in starbursts. 

Subject headings: stars: evolution — stars: atmospheres — cosmology: reionization 


1. INTRODUCTION 

Photons in the Lyman continuum (LyC) of hydrogen 
(energies E > 13.60 eV and wavelengths A < 911.753 A) 
play a crucial role in theoretical and observational studies 
of H II regions, active galactic nuclei (AGN), starburst 
galaxies, and ionization of the interstellar medium (ISM) 
and intergalactic medium (IGM). The dominant sources 
of ionizing photons in the ISM are 0-type stars (Spitzer 
1978), while the metagalactic ionizing background arises 
from both 0-type stars and AGN (Haardt & Madau 
2001, 2012; Shull et al. 1999). Photoionization by hot 
stars is also the likely process governing the H I epoch 
of reionization at redshifts z ss 7) ending the transition 
from the “dark ages” following cosmological recombina¬ 
tion at z ss 1100 through “first light” at z ss 30 — 50. 

The redshift frontier for detection of early galaxies is 
now at z ss 10 using the Lyman-break flux dropout tech¬ 
nique with the Wide-Field Camera-3 (WFC3) aboard the 
Hubble Space Telescope (Bouwens et al. 2011; Bradley 
et al. 2012; Ellis et al. 2013; Robertson et al. 2013; 
Schmidt et al. 2014). These distant galaxies allow us 
to explore stellar populations during the epoch (z = 
7 — 10) when reionization was well underway. The 
critical star formation rate (SFR) density, /sfr ~ 
0.02 Mq yr“^ Mpc“° at z « 7, defines the photoioniza¬ 
tion rate needed to balance IGM recombinations (Madau 
et al. 1999; Shull et al. 2012). However, estimating the 
number of ionizing photons in the radiation field requires 
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calibrating the 0-star production rates and extrapolating 
the galaxy luminosity function well below the observed 
populations (Trenti et al. 2011; Finkelstein et al. 2012; 
Robertson et al. 2013). The main parameters in psFR 
include the net production rate and escape fraction of 
Lyman continuum (LyC) photons by hot stars and struc¬ 
tural properties of the ISM and IGM. Gurrent estimates 
suggest that galaxies should be able to complete reion¬ 
ization by z ~ 7, although the details depend on LyC 
production efficiencies. 

The LyC production efficiency (ionizing photons per 
solar mass of star formation) is a standard parameter in 
cosmological simulations of galaxy formation and IGM 
reionization. Historically, many groups have assumed a 
constant rate of 10®°’^ photons s“^ per Mq yr“^. The 
current calculation explores the dependences of this pa¬ 
rameter on metallicity, initial mass function (IMF), and 
stellar rotation. Because of the importance of 0-star ion¬ 
izing photon production (Vacca et al. 1996; Schaerer & 
Vacca 1998; Oey et al. 2000; Kewley et al. 2001; Smith 
et al. 2002; Sternberg et al. 2003; Martins et al. 2005; 
Leitherer et al. 2010) to the study of starburst galax¬ 
ies and the epoch of reionization, there is a need for 
more accurate calculations. This can now be done, with 
improvements in stellar atmospheres modeling, includ¬ 
ing three-dimensional expanding atmospheres, non-LTE 
radiative transfer, populations of the hydrogen and he¬ 
lium ground states, and an increased number of spectral 
lines in the line blanketing and force modeling of stel¬ 
lar winds. For stellar populations, the LyC production 
rates depend on the initial mass function (IMF), evo¬ 
lutionary tracks, and model atmospheres, all of which 
may change at low metallicity and with rapid stellar ro- 
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tation (Schaerer 2003; Ekstrom et al. 2012; Georgy et al. 
2013). Some of the recent tracks also include the effects 
of rotation on stellar structure and evolution at different 
metallicities. 

In this paper, we compute the production efficiency 
(LyC photons per solar mass of star formation). We 
use new evolutionary tracks of various metallicities, with 
and without rotation, together with ionizing spectra from 
model atmospheres. We find the production rates of pho¬ 
tons in the ionizing continua of H I, He I, and He II 
by integrating over a grid of stellar tracks coupled to 
non-LTE, line-blanketed stellar atmospheres using the 
WM-basic code (Pauldrach et al. 2001; Puls et al. 2005) 
with a consistent metallicity. Section 2 describes our 
procedure for deriving the stellar LyC from evolutionary 
tracks and consistent model atmospheres. Section 3 ex¬ 
plains our methods for generating a grid of atmosphere 
models overlaid on evolutionary tracks of solar and sub¬ 
solar metallicities, with and without rotation. In Sec¬ 
tion 4 we present and analyze results for the ionizing 
photon luminosities and lifetime-integrated LyC produc¬ 
tion efficiencies, for both the hydrogen LyC as well as 
the continua of He I and He II (the latter is especially 
sensitive to metallicity and rotation). In Section 4.3, we 
apply the new LyC efficiencies to IGM reionization and 
critical SFR at z Ri 7. Section 5 provides a summary and 
a discussion of future directions in the field. 


2. COMPUTING THE STELLAR LYMAN CONTINUUM 


Our LyC calculation relies on coupling a full grid of 
stellar evolutionary tracks with model atmospheres. Al¬ 
though it would be desirable to explore a wider range 
of other parameters (mass-loss rates, clumpy and porous 
winds, and binaries) there are no complete evolution¬ 
ary tracks for such parameters. Thus, we rely on stellar 
grids from the Geneva group (Schaerer 2003; Ekstrom 
et al. 2012; Georgy et al. 2013). Our computational 
procedure begins with these evolutionary tracks, which 
provide the key stellar parameters: effective tempera¬ 
ture (Teff), surface gravity (logg), mass (M), and radius 
(i?*). From atomic data and stellar opacities, we run 
new model atmospheres to compute the spectral energy 
distribution (SED) with astrophysical flux {F\) in units 
of erg cm“^ s“^ A“^. We derive the integrated ionizing 
photon fluxes, qi (phot cm“^ s“^), from the star’s sur¬ 
face for the H I, He I, and He II continua by integrating 

the modeled spectra over wavelengths A < 


Qi = 
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Here, hc/X is photon energy and a|(,^ correspond to the 
H I, He I, and He II continuum edges at Aum = 911.753 A, 
504.259 A, and 227.838 A for indices * = 0 , 1 , and 2 
respectively. 

We combine the evolutionary tracks and grid of model 
atmospheres consistently to obtain an updated grid (Teg, 
log 5 ) where each point contains information on photon 
fluxes instead of the entire synthetic spectrum. The 
star’s ionizing photon luminosities Qi = AnR^qi (pho¬ 
tons s“^) are found from the star’s surface area at each 
point along the evolutionary track for a star of radius 
i?*(t) and initial mass M. The lifetime-integrated val¬ 


ues of Qo, Qi, and Q 2 are then derived for the H I, He I, 
and He II continua. Thus, (m) gives the total num¬ 
ber of ionizing (LyC) photons produced over a massive 
star’s life, where m = M/Mq is in solar-mass units. For 
0-type stars in the range 25 < m < 100, typical val¬ 
ues are ( 30 *°*^^ « (1 — 10) x 10®^ photons. These photon 
production numbers can then be summed over a stellar 
IMF using one of the common distributions (Salpeter 
1955; Kroupa 2001; Chabrier 2003) to determine the 
number of ionizing photons produced by a stellar pop¬ 
ulation. A similar grid of WM-basic model atmospheres 
at solar metallicity was computed by Sternberg et al. 
(2003). Our calculations are done with updated stellar 
parameters appropriate for the new evolutionary tracks. 

Reionization models and cosmological simulations re¬ 
quire a relationship between SFR and LyC production 
rate. This calibration is often expressed as the stan¬ 
dard ratio (Madau et al. 1999) of 10®®'^ photons s“^ 
per Mq yr“^. We adopt an average LyC production 
rate over the age of an OB starburst. Without great 
loss of accuracy, one can simply cancel the time units 
(3.1556 X 10^ s/yr) and quote the equivalent standard 
ratio as QlyC = 4 x 10®° LyC photons per solar mass of 
star formation. In this paper, we follow the convention 
used by Shull et al. (2012) in their reionization models, 
where they define the ratio 
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Here, ^(m) represents the differential distribution of stars 
per unit mass (not per logm as sometimes defined), 
and the integration is between mass limits mmin and 
Wmax- For convenient units, we express QlyC in 10®® 
photons per Mq of star formatior0. For a range of 
metallicities, 0.1 < (Z/Zq) < 1 , older model atmo¬ 
spheres, and a Salpeter IMF, Shull et al. (2012) found 
that QlyC ~ 0.0024— 0.0038 (xlO®® photons per Mq). 
Below, we update this calibration factor using new stel¬ 
lar data on LyC production and exploring the effects of 
lower metallicity and increased rotation. 


3. COMPUTATIONAL METHODS FOR TRACKS AND 
ATMOSPHERES 

The earliest stellar evolutionary tracks included only 
a small amount of stellar physics that was computation¬ 
ally allowed at the time. These models solved the equa¬ 
tions of stellar structure with little treatment of more ad¬ 
vanced physics. Some of the first standard evolutionary 
tracks came from the Geneva group in the early 1990’s 
(Schaller et al. 1992; Schaerer et al. 1993a,b; Charbonnel 
et al. 1993). After observations suggested higher mass 
loss rates for massive stars, additional set of tracks were 
produced (Meynet et al. 1994; Fagotto et al. 1994a,b). 
The next major advance was to add the effects of stel¬ 
lar rotation. Although the early efforts were productive, 
they only described some effects of rotation and did not 
extend to M < 9 Mq (Meynet & Maeder 2000). 


® The units of 10®® photons were chosen (Shull et al. 2012) to 
match the typical lifetime LyC photon production for O stars with 
mass 25 < m < 100. Integration over a stellar IMF reduces QLyC 
to much lower values, 0.002 — 0.004 since the low-mass stars con¬ 
tribute mass but essentially no LyC photons. 
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Fig. 1. — New rotating tracks (Ekstrom et al. 2012; Georgy et al. 
2013). Tracks computed at solar metallicity (top) and sub-solar 
metallicity (bottom). Note different post main sequence popu¬ 
lations at low metallicity, with fewer “blue-loop” excursions and 
fewer WR stars. 

For the past two decades, the most widely used stellar 
evolutionary tracks were those of Schaller et al. (1992). 
The most recent tracks (see Figure 1) computed at solar 
metallicity (Ekstrom et al. 2012) and sub-solar metallic¬ 
ity (Georgy et al. 2013) are expected to become the new 
standard in the field. These models surpass previous evo¬ 
lutionary tracks with their treatment of mass loss from 
the stellar wind, more accurate nuclear burning rates, 
and the addition of stellar rotation. Advances in helio¬ 
seismology and other stellar observations have improved 
our understanding of the interior properties of rotating 
stars. From these data, Ekstrom et al. (2012) created 
evolutionary tracks and isochrones for a wide mass range 
of rotating and corresponding non-rotating stars. These 
models can be used to create an initial picture of rotating 
stellar evolution. One of the most important differences 
with rotation is mixing caused by meridional circulation, 
which continuously provides new fuel to burn in the core 
and increases the lifetime of a star. 

Improvements in the quality of model atmospheres in¬ 
clude non-LTE treatment of the radiation field, wind 
and line blanketing, and increased code efficiency. The 
first upgrades involved an efficient method for approx¬ 


imating the radiation pressure due to spectral lines. 

The CAK model (Castor, Abbott & Klein 1975) aver¬ 
aged the radiative force over all lines instead of treat¬ 
ing them separately. This increase in radiative force 
led to higher mass loss rates. Even with this advance, 
atmosphere models were restricted to local thermody¬ 
namic equilibrium (LTE), requiring a blackbody radia¬ 
tion field and Boltzmann level populations. Pauldrach 
et al. (1986) upgraded the CAK model by neglecting the 
radial-streaming approximation for photons that drive 
the stellar wind. Rotational effects were found to be 
negligible for rotational velocities Urot < 200 km s“^. 

Advances in model atmospheres have been more con¬ 
sistent than advances in evolutionary tracks, primarily 
because of the availability of fast computers. The most 
efficient codes still use a modified version of the CAK 
method to compute line forces, because it is fast and 
its approximations provide efficiency. Perhaps the most 
important improvement is the extension of the models 
into the non-LTE regime, allowing the radiation field 
to deviate from Planckian and level populations from 
Boltzmann. Other advances include solving the radia¬ 
tive transfer equations in co-moving coordinates, mov¬ 
ing away from plane-parallel models, and increasing the 
accuracy of atomic data. With these recent improve¬ 
ments, a model atmosphere can now be computed in less 
than ten minutes, allowing thousands of models to be 
computed over the course of the project and providing 
accurate time resolution for the evolution of stellar spec¬ 
tra. A few of the most prominent atmosphere codes 
are: TLUSTY (Lanz & Hubeny 2003, 2007), FASTWIND 
(Puls et al. 2005), CMFGEN (Hillier & Miller 1998), and 
WM-basic (Pauldrach et al. 2001). 

For our model atmospheres we chose the WM-basic 
cod43 primarily because of its hydrodynamic solution of 
spherically expanding atmospheres with line blanketing 
and non-LTE radiative transfer and its treatment of the 
continuum and lines in the ionizing EUV. This code is 
also a standard for theoretical atmosphere modeling in 
the population synthesis code Starburst99 (Leitherer 
et al. 1999, 2014). An extensive discussion of hot-star 
atmosphere codes appears in papers by Martins et al. 

(2005) and Leitherer et al. (2014). A comparison of 
the ability of these codes to fit stellar lines, effective 
temperatures, and surface gravities was made by Puls 
(2008) and Massey et al. (2013). In Figure 2 we show 
the results of WM-basic and CMFGEN for three cases with 
Teff ~ 30kK, 40kK, and 50kK. The agreement is reason¬ 
ably good, given the different model input parameters 
required to arrive at similar effective temperatures, but 
slightly different logg values. 

3.1. Evolutionary Tracks 

We use evolutionary tracks from Schaller et al. (1992), 

Ekstrom et al. (2012), and Georgy et al. (2013). These 
models include fixed prescriptions for stellar mass loss 
during late stages of evolution, and they provide a con¬ 
venient means of assessing the effects of metallicity and 
rotation. It would be useful to assess the dependence 
of LyC production rates on alternative prescriptions for 
mass loss, particularly during late stages of massive- 
star evolution. However, at the current time, we are 

^ Code can be found at http://www.usni.uni-muenchen.de/people/adi/Prograi 
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Fig. 2.— Overlaid comparison of model atmospheres at effective temperatures ^ 30 kK, 40 kK, and 50 kK, computed using the 
WM-basic and CMFGEN code with emergent fluxes F\ in units 10^^ erg cm“^ s“^ A“^. Because of slight differences in mass-loss rates and 
required input parameters for the stellar models, we are unable to precisely match effective temperatures (Teg) and surface gravities (logg). 
The WM-basic model is shown in black and the CMFGEN model in red. 


































































































































5 


constrained to using the available evolutionary tracks. 
The field of massive-star evolution continues to present 
challenges to the correct implementation of these effects 
(Vink 2014; Martins 2014; Hirschi 2014) as well as those 
of binaries and stellar mergers (Eldridge et al. 2008; 
de Mink et al. 2013, 2014). Populations of very mas¬ 
sive stars (150-300 Mq) may exist in the Arches cluster 
(Figer et al. 2002; Crowther et al. 2010), in the dense 
R136 star cluster in 30 Doradus, and in the Galactic 
cluster NGC 3603 (Growther et al. 2010; Martins 2014). 
Figer (2005) used the Arches cluster to set an upper mass 
limit of 150 Mq, owing to an absence of star^l with ini¬ 
tial masses > 130 Mq. These mass estimates are subject 
to uncertainties in stellar luminosity arising from cluster 
distance, extinction laws, and effective temperatures. As 
summarized in recent reviews (Vink 2014; Martins 2014), 
there are several stars that could have masses between 
100-200 Mq. Such stars could provide additional sources 
of LyG by extending the IMF above the last grid point 
on the evolutionary tracks. On the other hand, the most 
massive stars in dense clusters are deeply embedded in 
star-forming clouds, and their LyC may not escape until 
after these regions clear of gas and dust. 

For the current calculation, we truncate all IMFs at 
a maximum mass of 120 Mq, and we use the LyC pro¬ 
duction rates to describe massive star formation regions 
and reionization of the IGM. The tracks are computed 
at five different metallicities, of which we consider three 
for non-rotating stars. In the Hertzsprung-Russell (H- 
R) diagram stellar evolution defines the basic parameters 
{M*, L*, Teff, \ogg) which are related by self-consistency. 
For example, surface gravity g = GM^/Rl depends on 
stellar mass M* and radius i?*, and stellar luminosity 
L* = depends on effective temperature T^s 

and i?*. The ionizing fluxes are sensitive primarily to 
Teff. We use different evolutionary tracks to investigate 
how different stellar physics can affect the production of 
LyC photons. The first set of tracks are from Schaller 
et al. (1992), a past standard for non-rotating stars with 
sub-solar, solar, and super-solar metallicities, Z = 0.004, 
0.02 and 0.04, respectively, by mass. The second set 
of tracks comes from the same group (Fkstrbm et al. 
2012) with a mass range extended to 0.8 Mq < m < 
120 Mq and finer mass resolution. One essential ad¬ 
vance is their treatment of rotation. We use two sets of 
tracks with metallicities Z = 0.004 and Z = 0.014, each 
computed with an initial rotation of 40% of the critical 
rotation speed. We also consider their tracks without 
rotation. Since the Schaller et al. (1992) work, the ac¬ 
cepted value for the solar metallicity has decreased to 
either Zq = 0.014 (Asplund et al. 2009) or Zq = 0.0153 
(Caffau et al. 2011). Fkstrom et al. (2012) calculated 
both rotating and non-rotating tracks at solar metallicity 
{Z = 0.014), while Georgy et al. (2013) computed both 
rotating and non-rotating tracks at sub-solar metallicity 
{Z = 0.002). The Geneva group has not yet provided 

® Subsequent studies (Espinoza et al. 2009) of stars in the Arches 
cluster found a classical IMF with mass distribution dN/d(log m) oc 
yyj-l.l± 0.2 (gaipeter index is F = 1.35) for masses m > 10 Mq 
and measured through nearly 20 magnitudes of visible extinction. 
Crowther et al. (2010) examined several of the Arches cluster stars 
and derived larger luminosities owing to a larger assumed distance 
and different extinction law. They reported a few masses in excess 
of 160 Mq. 


new tracks with super-solar metallicity. 

3.2. Model Atmospheres 

After compiling the evolutionary tracks, we produce 
model spectra at each time step along the track. We use 
the atmosphere code WM-basic because of its efficiency 
and ability to generate accurate spectra, including lines 
in the UV and FUV. The code requires only a few ba¬ 
sic parameters (Tes, logg, i?*, abundances) to compute a 
model. For each set of tracks we compute a grid of atmo¬ 
spheres in the (Tea, log g) plane and follow the evolution 
of the ionizing spectra of a star through time. We use the 
flux distributions to find the total LyG flux and integrate 
to calculate the ionizing photon production for each ini¬ 
tial mass. An atmosphere model can be completed in 
10-15 minutes on a desktop computei0. Over the course 
of our project, we ran hundreds of atmosphere models 
to achieve accurate time resolution along evolutionary 
tracks. One full grid of models for a set of evolutionary 
tracks contains around 100 atmosphere models. 

The process of model atmosphere calculation comes 
in three different pieces: the hydrodynamics, calcula¬ 
tion of the occupation numbers and radiation held, and 
production of the synthetic spectrum. The hydrody¬ 
namics section is where the initial parameters are input, 
most importantly Tea, logg, and i?*. From these, the 
code hnds the mass, luminosity, and other basic data. 
The next input are the abundances of elements between 
H and Zn, with three options: solar metallicity, solar 
metallicity while dehning the H and He abundances, and 
custom-dehned abundances. Because the default “solar” 
abundances programmed into the code were outdated, 
we used the recent standard of Zq = Q .014 by mass (As¬ 
plund et al. 2009). The code requires input of the metal 
abundances relative to solar, which we re-examine after 
a model is run. In its normalization procedure, one can 
alter the abundances. The resulting values may differ 
from the input values, but generally by only a few per¬ 
cent. Once the input file is written, the model can be 
run using the graphical user interface. 

For the radiation force on the atmosphere and wind, 
we use three line-force parameters {k, a, 6) in the GAK 
method. This method approximates the sum of a large 
number of spectral lines, with temperature assumed to 
be constant throughout the atmosphere to increase the 
efficiency. The hydrodynamics calculation is iterated 
two more times. The hrst time, the temperature gra¬ 
dient is computed using the continuum opacities. In 
the last iteration the temperature is computed including 
the line opacities. The code then computes the photon 
energy, occupation numbers, and radiation held in the 
NLTF regime. This step provides the radiation held, h- 
nal temperature structure, information about opacities, 
and occupation numbers for the various elements us¬ 
ing detailed atomic models. The hnal step in the cal¬ 
culation provides the synthetic spectrum. Figure 3 il¬ 
lustrates the continuum shapes from models with Teff 
ranging from 31,800 K to 72,300 K, particularly the con¬ 
tributions of emission lines, the depth of the Lyman 
edge, F(912-)/F(912+), and the ratio, F(1500)/F(910), 

® We used an AMD Phenom II X6 1075T 3.0GHz processor and 
32Gb RAM. With the multiple processor cores of this machine, we 
typically run up to six models at a time. 


erg s-‘cm-U“‘ (10“ ) Fj erg s ‘cm "A ‘ (10“) Fj ergs-'cm-^ ‘ (10“) 
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Fig. 3.— Six model atmospheres computed using the WM-basic atmosphere code and taken from the Ekstrom et al. (2012) grid of rotating 
tracks at solar metallicity (Z = 0.014). Models are shown in lOkK intervals in Tefr; we include 45kK because of the rapid change in strength 
of the H and He edges in that temperature range. Models are all for 0-type stars near the main sequence. 
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of fluxes at 1500 A and 900 A. The latter two ratios 
have been used to calibrate the escape fraction (/esc) 
of LyC radiation. This observational technique was in¬ 
troduced by Steidel et al. (2001) and employed by In- 
oue et al. (2005), Shapley et al. (2006), and Mostardi 
et al. (2013) to estimate /esc in star-forming galaxies at 
0 = 2.85 — 3.27. In their method, the observed flux- 
density ratio, [F(1500)/T'(910)]obs, is compared to esti¬ 
mates of the intrinsic ratios determined from the models. 
We return to this issue in Section 4.5, where we dis¬ 
cuss the intrinsic values, [T'(1500)/F(910)]int ~ 0.4 —0.7 
found in our modeled composite spectra of OB associa¬ 
tions with coeval starbursts. 

In the models of Figure 3 with Tes < 50,000 K, the 
absorption edge at the Lyman limit is clearly visible at 
A = 912 A. Figure 4 shows two models at Teg w 41, 000 K 
with metallicities Z = 0.014 (solar) and 0.004 (sub-solar) 
by mass. Each has an intrinsic Lyman edge resulting 
from the non-LTE level population of the ground state 
of H I. Differences in metallicity have a minor effect on 
the depth of the edge, which arises primarily from sur¬ 
face abundances and heating. However, the EUV fluxes 
are affected by emission lines in the stellar wind, which 
increase in strength at higher metallicity (left panel of 
Figure 4). 

3.3. Grid of Atmosphere Models for Evolutionary 
Tracks 

In order to find the integrated LyC photon production 
over an 0-star’s lifetime, we set up a grid of spectra in the 
(Teff,log( 7 ) plane defined by the initial mass of the evo¬ 
lutionary track (Figure 5). We distribute the grid points 
evenly along points of equal time. Because the grid has 
a higher density where the stars spend more time, we 
obtain a higher accuracy of photon fluxes over the stellar 
lifetime and improve the accuracy of the integrated pho¬ 
ton fluxes. As the most massive stars shed their atmo¬ 
spheres in late stages of evolution (Figure 6), they move 
quickly through “blue loops” and into the Wolf-Rayet 
(WR) phase. Although the most massive stars spend 
only a short time in these regions of the H-R diagram, 
their LyC production rate is high. Therefore, we com¬ 
pute extra atmosphere models and ensure that we place 
sufficient grid points during these late stages to capture 
the ionizing radiation. 

The WR stars deserve further discussion as the evolu¬ 
tionary descendants of massive stars (M > 25 Mq) char¬ 
acterized observationally by broad emission lines, strong 
stellar winds, and mass-loss rates up to I0“® Mq yr“^ 
(Maeder & Conti 1994; Crowther 2007; Massey 2013). 
Both luminous {L > 10® Lq) and hot (25,000- 

100,000 K), WR stars can contribute substantially to the 
LyC production of a stellar population, particularly in 
the He I and He II continua. However, the canonical esti¬ 
mate that WR stars are 10% of the massive-star popula¬ 
tion is based on little more than the He-burning lifetimes 
of stars with M > 25 Mq estimated from uncertain evo¬ 
lutionary tracks. Because 0-star mass-loss rates depend 
on metallicity, it should be easier to become a WR star at 
high metallicity. This expectation appears to agree with 
observations, which show fewer WRs at low metallicity in 
Local Group galaxies (Maeder et al. 1980; Hamann et al. 
2006; Neugent et al. 2012; Massey et al. 2014). However, 


the progenitor masses and details of post-main-sequence 
stellar evolution that lead to the various (WC, WN, WO) 
sub-classes are still poorly understood. Recent surveys 
of WN and WC stars in the Galaxy, LMC, and M31 
(Sander et al. 2012; Sander et al. 2014; Hainich et al. 
2014) suggest that WN progenitors originate from initial 
masses between 20-60 Mq. The ratios of WR stars to 
O stars (and of WN and WC subtypes) are frequently 
used as tests of massive-star evolutionary models. Ow¬ 
ing to uncertainties in distance and reddening and se¬ 
lection biases in detecting WR stars relative to 0-type 
stars, astronomers probably cannot measure the Galactic 
(WR/0) ratio to better than 50% (Massey et al. 1995). 

Current surveys (Mauerhan et al. 2011; Shara et al. 
2012) identified ~ 500 WR stars within 3-4 kpc of the 
Sun (plus the Galactic center). From Galactic distri¬ 
bution models, Shara et al. (2009) estimated that the 
Milky Way could have up to 6500 WR stars. A more re¬ 
cent study (Rosslove & Growther 2015) has revised this 
estimate downward to 1900 ± 250. Maeder & Meynet 
(2005) found WR/0 « 0.02 in models with no stellar 
rotation, increasing to WR/0 ~ 0.05 — 0.07 for rotation 
speeds of 300 km s“^ and metallicities Z = 0.01 — 0.02. 
The WR/0 ratio drops to only a few percent in the 
Magellanic Clouds (Meynet & Maeder 2005), although 
the observed WR/0 ratios continue to change with new 
LMC/SMC surveys (Massey et al. 2014). Georgy et al. 
(2012) compared recent observations of massive-star pop¬ 
ulations WR/0 « 0.12 ± 0.03 in the Milky Way with 
their solar-metallicity models {Z = 0.014) with and with¬ 
out rotation. They found that models with rotation 
could generate WR/0 w 0.07, and they suggested that 
close-binary interactions might explain the remainder. 
Thus, we believe that our calculations provide a con¬ 
servative estimate of the LyG produced by O stars and 
their post-main-sequence successors, consistent with the 
latest Geneva tracks. At the low metallicities expected in 
the high-redshift galaxies responsible for reionization, the 
WR populations will likely be much smaller. Future cal¬ 
culations may improve the accuracy of the populations 
by including effects of binaries and improved prescrip¬ 
tions for stellar mass loss. 

Our overall goal is to compute the ionizing photon pro¬ 
duction value at each of the timesteps along the tracks. 
Because running a model at every time step would re¬ 
quire thousands of runs, we chose a different approach. 
We ran model atmospheres along ~ 100 key points and 
then interpolate to define the grid of spectra for each 
track. Typically, we have a few tens of models for each 
initial mass, totaling ~100 models for each set of tracks. 
For blue loops in two most massive stars (85 Mq and 120 
Mq) the stars move far to the left in the H-R diagram, 
Teff > 70 kK. This is too hot for our atmosphere code to 
achieve convergence, as the star exceeds the Eddington 
luminosity (with full opacity sources). Instead of interpo¬ 
lating, we match the point along the evolutionary track 
with the closest model atmosphere. The error of this ap¬ 
proximation is small because of the short time the stars 
spend in these regions. For a general point on a track 
we choose four model atmospheres on the grid surround¬ 
ing it and linearly interpolate to the track across the 
parallelogram connecting them in (Tgff, log) space. Two 
of the four points are along the same track and provide 
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Fig. 4.— Two model atmospheres with Teff ~ 41kK, computed using the WM-basic atmosphere code showing metallicity effects on the 
strength of the Lyman edge at 912 A and emission lines from the stellar wind. Left: Solar metallicity model {Z = 0.014) has a 57% 
drop in flux at the continuum edge and substantial metal emission lines in the EUV band (300-900 A). Right: Sub-solar metallicity model 
(Z = 0.004) has a 68% continuum drop and many fewer emission lines. 




Fig. 5.— Top: Evolutionary tracks (Ekstrom et al. 2012) at solar 
metallicity computed with rotation and labeled by initial mass (in 
Mq). Bottom: Corresponding non-rotating tracks. Points along 
tracks show values of Teff and L at which model atmospheres were 
computed. 



Fig. 6. — The LyC photon luminosities, Qo{t), versus time for 
Ekstrom et al. (2012) rotating evolutionary tracks with Z = 0.014. 
Each line shows a track with different initial mass from 120 Mq 
to 25 Mq (top to bottom). Rapid spikes show “blue loops” and 
Wolf-Rayet phase at the end of the star’s lifetime. Although WR 
stars take up a small fraction of the lifetime, we have ensured 
that sufficient model atmospheres are placed to capture the LyC 
produced in these late phases of evolution. 

the main part of the interpolation. The two points out¬ 
side the track provide minor corrections. Some regions 
of the grid that do not allow us to use all four points for 
interpolation. In these regions, we use a smaller num¬ 
ber of suitable points to complete the process. There are 
also times when the evolutionary tracks move outside the 
grid. We then find the closest grid point and adopt its 
value for the photon production. We then integrate the 
ionizing photon luminosities for each initial mass from 
the evolutionary tracks to find the total number of ion¬ 
izing photons produced over the stellar lifetime. 
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Fig. 7. — Lifetime-integrated number of H-ionizing photons, 
(m), for three older evolutionary tracks, labeled “Schaerer” 
and taken from Schaller et al. (1992) at three different metallici- 
ties. For comparison, we show LyC production for two new solar- 
metallicity tracks (Ekstrom et al. 2012) computed with and without 
rotation. 



Initial Mass (M^) 

Fig. 8. — Comparison of lifetime-integrated LyC photon produc¬ 
tion, Qg*°*^(m), for new evolutionary tracks (Ekstrom et al. 2012; 
Georgy et al. 2013), with and without rotation, at solar (Z = 0.014) 
and sub-solar (Z = 0.002) metallicity. 

4. RESULTS 

4.1. Lifetime-Integrated HI Photon Luminosities 

Using our coupled tracks and model atmospheres, 
we derive the total hydrogen LyC photon production, 
over the lifetime of individual massive stars 
(Table 1). Figure 7 shows these values for three older 
evolutionary tracks (Schaller et al. 1992) with no rota¬ 
tion and at three different metallicities {Z = 0.004, 0.02, 
0.04). These are compared to two new tracks (Ekstrom 
et al. 2012; Georgy et al. 2013) at solar metallicity, Z = 
0.014, with and without rotation. Clearly, the highest 
LyC production comes from the rotating tracks, which 
are typically 50% higher, with = (3 — 11) x 10®^ 

photons for initial masses of 40-120 Mq. Models calcu¬ 
lated at lower metallicity consistently produce more ion¬ 
izing photons, although the increase is smaller (^ 10%) 
than in rotating models. 


Figure 8 compares for the new rotating and 

non-rotating tracks (Ekstrom et al. 2012; Georgy et al. 
2013) at the two available metallicities, Z = 0.014 and 
Z = 0.002. Here again, we see that LyC production is 
elevated by ^ 50% for models with rotation over most 
of the mass range of 0-stars. Our integrations with the 
new non-rotating tracks produce QQ°^\m) curves (two 
dashed lines in Figure 8 ) that are nearly linear with mass, 
while production rates with the older tracks (Schaller 
et al. 1992) flatten out for the 85 Mq and 120 Mq 
models. Interestingly, the new rotating-track models of 
Georgy et al. (2013) show considerably more LyG pho¬ 
tons at 120 Mq, continuing the linear trend of QQ°^\m). 
However, those tracks were computed at low metallicity 
{Z = 0.002) which may explain part of the increase in 
LyC production. 

In general, stars with lower metallicity have a reduced 
rate of CNO-cycle burning, which results in a contrac¬ 
tion and heating of the star’s core and increased burning 
through the proton-proton chain. Consequently, the sur¬ 
face temperature is hotter and the star produces more 
ionizing photons. Additionally, fewer metals in the at¬ 
mosphere produce less line blanketing, which alters the 
spectral shape in the UV and EUV. The effects of stellar 

rotation on LyC production and QQ°*\m) arise partly 
from the prolongation of hydrogen burning, as new fuel 
is mixed into the core. Rotating models are more lu¬ 
minous, because of their larger convective cores, and 
their surfaces are hotter owing to enhanced helium abun¬ 
dances. As our calculations bear out, the ionizing fluxes 
are greater, as are their lifetime integrated LyC pro¬ 
ductions. Thus, owing to the complicated evolutionary 
tracks of massive stars through different ranges of pa¬ 
rameters (Teff, log 5 , L, i?*) there are multiple factors 
contributing to the increase of LyC production found at 
low metallicity and with enhanced rotation. 

4.2. Integration Over Initial Mass Functions 

With these data and the IMF of a stellar population, 
we can calculate the number of LyC photons produced 
per solar mass of star formation. To do this, we in¬ 
tegrate the IMF, multiplied by the lifetime inte¬ 

grated photon number Qi{m) for the three ionization 
continua {i = 0,1,2). We also consider three differ¬ 
ent IMFs. The hrst is adapted from Salpeter (1955), 
who found ^(m) = with index a = 2.35 for 

stars between 0.4 Mq and 10 Mq. We extend the mass 
range to 0.1 < m < 120 (to in Mq) and take K as 
a normalization constant. The Salpeter function is the 
steepest of the three IMFs, maintaining the same in¬ 
dex through its entire mass range, whereas recent sur¬ 
veys suggest a turnover at low mass. The second IMF 
is from Kroupa (2001), a piecewise-continuous power 
law, split into three mass ranges, each with a differ¬ 
ent value of a. We integrate over the two stellar mass 
ranges between 0.08 Mq and 120 Mq, in which a = 1.3 
(0.08 < TO < 0.5) and a = 2.3 (0.5 < to < 120). The 
third IMF is from Chabrier (2003), with the same index 
{a = 2.3) for stars with to > 0.5. For 0.08 < to < 0.5 
this IMF has the more complicated log-normal form; 
^(to) = Aexp[—(logTO — logTOc)^/2cr^]. 

Table 2 summarizes the results of our calculations for 
various combinations of IMF, metallicities, and evolu- 
















10 


tionary tracks, with and without rotation. The LyC 
production efficiency is defined as the lifetime-integrated 
number of ionizing photons per solar mass of star forma¬ 
tion. We list three versions of these parameters, QhyC, 
Quei and QHeii, corresponding to the ionizing continua 
of H I, He I, and He II. These are sensitive to the sets 
of model atmospheres and stellar tracks, as well as the 
IMF index (a) at the high-mass end, and the integrated 
mass range (wniin and Wmax)- 


4.3. Reionization and LyC Production Effiicency 

We now apply the new stellar models to reionization by 
comparing our calculated values of QhyC to the LyC pro¬ 
duction efficiencies calculated by Shull et al. (2012). The 
ionizing photons QhyC produced per Mq of star forma¬ 
tion are important for determining when IGM reioniza¬ 
tion was mostly complete. This can be assessed by find¬ 
ing the critical star formation rate density (Madau et al. 
1999) necessary to keep the IGM ionized, balancing pho¬ 
toionizations with hydrogen recombinations. The critical 
density was recently evaluated (Shull et al. 2012) using 
updated estimates for LyG production rates (Quyc) and 
LyC escape fraction (/esc) as well as IGM structure (den¬ 
sity and temperature). Hydrogen (case-B) recombina¬ 
tion rates depend on electron temperature as and 

are enhanced by a “clumping factor” Ch = 

With calculations of IGM clumping, escape fraction, and 
LyC production efficiencies, Shull et al. (2012) suggested 
fiducial values of Ch ~ 3, /esc ~ 0.2, Te = (10^ K)r 4 , 
and LyC efficiencies QhyC = 0.004 (in 10®® photons per 
Mq) to find 


PSFR = (0.018 Mq yr ® Mpe ®) 
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gH/3 \ / 0.004 \ 

\QLycJ " 


(3) 


The value of /esc for low-mass galaxies with a cloudy ISM 
has been estimated to be 10-20% (Fernandez & Shull 
2011). Depending on choices of Ch and T^, the escape 
fraction needs to be at least 10% in order to complete 
IGM reionization by z = 7 (Shull et al. 2012; Robertson 
et al. 2013; Finkelstein et al. 2013). 

In our study, we have explored three factors that af¬ 
fect the LyG production efficiency: (1) decreases in stel¬ 
lar metallicity; (2) increases in stellar rotation; (3) shape 
of the IMF. From our new calculations of ionizing pho¬ 
ton production rates we assess how LyC production ef¬ 
ficiency is affected by choices of stellar tracks (rotation, 
metallicity), model atmospheres, and IMF. By compar¬ 
ing results for various combinations of these parameters, 
we findthe general trends with metallicity and rotation. 
We currently have the full range of evolutionary tracks 
at two metallicities, solar and sub-solar. We integrate 
over three IMFs, labeled in Table 2 as (1,2,3) corre¬ 
sponding to (Salpeter, Kroupa, Chabrier). The stellar 
tracks are labeled as (4,5,6) corresponding to (Schaller 
1992, Ekstrom 2012, Georgy 2013). The latter two tracks 
come with and without rotation. Our results (Figures 
7 and 8) show that for each set of non-rotating evolu¬ 
tionary tracks their rotating counterparts produce more 
photons. Metallicity is less important for the overall ef¬ 
ficiency, QhyC, although it can be a factor at specific 


masses. For the 40 Mq and 60 Mq models with rotation, 
the solar metallicity models produce more photons than 
the sub-solar models. Owing to the power-law decrease 
of the IMF at high masses, this small difference in ion¬ 
izing photon production results in the solar-metallicity 
models producing slightly more photons per solar mass. 
For fixed evolutionary tracks, metallicities, and model 
atmospheres, adopting the Kroupa or Chabrier IMF in¬ 
creases the LyC efficiency by factors of 1.56 and 1.64, 
respectively, relative to the Salpeter IMF. 

We now perform a detailed comparison of our com¬ 
puted values of QhyC to the fiducial value QhyC = 0.004 
suggested by Shull et al. (2012). With different combi¬ 
nations of IMFs and evolutionary tracks, QhyC extends 
from 0.00315 to 0.00940 for solar metallicity and from 
0.00402 to 0.00910 for sub-solar metallicity (Table 2). 
This spread in efficiencies translates to a factor of 3 in 
critical SFR density. Values of QhyC computed with ro¬ 
tating evolutionary tracks and the Kronpa IMF were at 
the upper end of the range. One can further analyze 
these Q-values to separate the effects of reduced metal¬ 
licity from those of rotation. The three IMFs only change 
the overall normalization, because of their different mass¬ 
loading by low-mass stars. Because the Salpeter IMF in¬ 
cludes more low-mass stars, the absolute values of QhyC 
are reduced; these stars add mass but produce few ion¬ 
izing photons. 

First, we examine the effects of metallicity. For the 
“old” non-rotating tracks (Schaller et al. 1992) at solar 
metallicity {Z = 0.020) we find QhyC = 0.00367, 0.00572, 
and 0.00370 for the three IMFs (1,2,3) considered here 
(Salpeter, Kroupa, and Ghabrier, respectively). For the 
new non-rotating tracks (Ekstrom et al. 2012) at sub¬ 
solar metallicity {Z = 0.004) we find QhyC = 0.00399, 
0.00623, and 0.00402 for IMFs (1,2,3). Taking the ap¬ 
propriate ratios for the old tracks and corresponding 
model atmospheres, we find that decreasing the metal¬ 
licity from solar (Z = 0.02) to sub-solar {Z = 0.004) 
increases the LyC efficiencies by only 9%. That is, the 
ratio QuyC (subsolar)/QuyC (solar) = 1.09 for all three 
IMFs. However, performing the same calculation for the 
new non-rotating tracks (Ekstrom et al. 2012; Georgy 
et al. 2013), we find a 34% increase at sub-solar metal¬ 
licity, with QuyC (subsolar)/QuyC (solar) = 1.34. For 
the new rotating tracks, the LyC efficiencies at solar 
and sub-solar metallicity are essentially the same, with 
QuyC (subsolar)/QLyC(solar) = 0.97 for all three IMFs. 
Chabrier and Kroupa IMFs have the same power-law 
slope for high mass (a = 2.3), but each uses a differ¬ 
ent fit for low-mass stars. The Chabrier IMF has more 
low-mass stars compared to the Kroupa IMF and there¬ 
fore more stellar mass that is not producing ionizing 
photons. Next, we examine the effect of rotation on 
QhyC at fixed metallicity. For all three IMFs, the evolu¬ 
tionary tracks with rotation produce more ionizing pho¬ 
tons per solar mass. At solar metallicity {Z = 0.014) 
the Ekstrom et al. (2012) tracks and corresponding 
model atmospheres yield 87% higher LyC efficiency, with 
QuyC (i'ot)/(5LyC (nonrot) = 1.87. At sub-solar metallic¬ 
ity {Z = 0.002) the Georgy et al. (2013) tracks and cor¬ 
responding atmospheres were 35% more efficient in LyC 
production, with QLyc(i'ot)/(3LyC(nonrot) = 1.35. 

As a general conclusion, we find that increased stellar 
rotation is more important than reduced metallicity in 
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boosting the LyC efficiency. By comparing subsets of 
tracks, metallicities, and IMFs, we can recover more spe¬ 
cific results. We label the LyC production efficiencies by 
QuyC (IMF, Tracks), where IMFs and tracks are defined 
as in Table 2 (1 = Salpeter, 2 = Kroupa, 3 = Chabrier) 
with non-rotating tracks labeled by (4 = Schaller; 5 = 
Ekstrom) and rotating tracks by (5 = Ekstrom, 6 = 
Georgy). Eor the relevant combinations of parameters, 
we find the following efficiencies for “non-rotating, solar 
metallicity” stellar populations: 

QLyc(l,4) = 0.00367 QLyc(2,4) = 0.00572 
QLyc(3,4) = 0.00370 (Z = 0.020) 
QLyc(l, 5) = 0.00322 QLyc(2,5) = 0.00503 

QLyc(3,5) = 0.00325 {Z = 0.014) (4) 

The range of QhyC between 0.00322 — 0.00572 arises pri¬ 
marily from the choice of IMF, which affects the “mass¬ 
loading” at the low end. Next, we make a similar tabula¬ 
tion for “non-rotating, sub-solar metallicity” populations 
(either Z = 0.004 or 0.002): 

QLyc(l,4) = 0.00399 gLyc(2,4) = 0.00623 
QLyc(3,4) = 0.00402 (Z = 0.004) 
QLyc(l,6) = 0.00433 QLyc(2,6) = 0.00676 

QLyc(3,6) = 0.00436 {Z = 0.002) (5) 

As before, the range of QhyC from 0.00399 — 0.00676 
arises from the assumed IMF. Finally, for the rotating 
tracks at both solar and sub-solar metallicities, we find: 


QLyc(l, 5) = 0.00603 QLyc(2,5) = 0.00940 
QLyc(3,5) = 0.00607 (Z = 0.014) 
QLyc(l, 6) = 0.00583 gLyc(2,6) = 0.00910 

QLyc(3,6) = 0.00588 (Z = 0.002) (6) 


The range of QhyC from 0.00583 — 0.00940 arises mainly 
from the assumed IMF. 

Adopting an IMF with fewer low-mass stars or us¬ 
ing tracks with increased stellar rotation produces larger 
boosts in LyC efficiency than sub-solar metallicity. For 
non-rotating stars at solar metallicity, a typical value is 
QhyC = 0.004 ± 0.002, as found in our previous study 
(Shull et al. 2012). Including stellar rotation could in¬ 
crease the efficiency to QhyC = 0.0075 ±0.0015, although 
those stellar tracks assumed a single value of rotation at 
40% of breakup, which may be excessively high. Allow¬ 
ing for the range in QhyC efficiencies produced by as¬ 
sumptions about IMF, metallicity, and rotation, we rec¬ 
ommend adopting an overall LyC production efficiency 
of QhyC = 0.006 ± 0.002 in Equation (3), or a photon 
production calibration of (6 ± 2) x 10®° LyC photons per 
Mq of star formation. This is equivalent to a SFR cal¬ 
ibration of photons s“^ for SFR in Mq yr“^ 

and represents a 50% increase over previous estimates. 
With the new LyC rates, the critical SFR to maintain 
reionization can be rewritten 


PSFR = (0.012 Mq yr ^ Mpc ®) 
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Fig. 9.— Reionization histories and critical SFRs for two LyC 
production efficiencies, computed with our reionization simulator 
(Shull et al. 2012). Left: Ionized volume fraction Qi and inte¬ 
grated CMB optical depth, Te, to electron scattering for two LyC 
production rates: QLyC = 0.006 (solid red curve, model A) and 
QLyC = 0.004 (dashed blue curve, model B), both in 10®^ pho¬ 
tons per Mq of star formation. The dot-dashed curve at 2 < 5 
follows the He II history, governed by quasars and completed at 
z 3. Right: Critical SFR densities required to maintain IGM 
ionization compared to observed and extrapolated SFR (green 
dot-dashed curve). Reionization models (assume clumping factor 
Ch = 3, LyC escape fraction /esc = 0.2, and electron temperature 
Te = 20, 000 K. The SFR history comes from the galaxy luminosity 
function (Trenti et al. 2010) extrapolated down to absolute magni¬ 
tude Mab = —10 in rest-frame UV and converted to SFR density 
using the prescription in Madau et al. (1999). 


The 50% boost in LyC production efficiency is an im¬ 
portant ingredient in assessing whether IGM reionization 
was complete by z w 7, as suggested by recent observa¬ 
tions (Oesch et al. 2014; McLure et al. 2013; Robertson 
et al. 2013; Finkelstein et al. 2012). 

These considerations are illustrated in Figure 9, which 
shows the reionization histories and critical SFR den¬ 
sities derived with the reionization simulator developed 
in a previous study (Shull et al. 2012). The left panel 
shows the redshift evolution of the ionized volume frac¬ 
tion, Qi{z), for two values of LyC production efficiency, 
QhyC = 0.004 (the old standard) and QuyC = 0.006 (the 
new efficiency). Accompanying curves show the evolu¬ 
tion of CMB optical depth, Te, to electron scattering, an 
integral constraint on the ionization history. These two 
models, labeled A and B, adopt identical physical pa- 
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rameters, Ch = 3, /esc = 0.2, and = 20, 000 K. The 
SFR history follows that of Trent! et al. (2010) with the 
galaxy luminosity function extrapolated down to abso¬ 
lute magnitude Mab = —10 in the rest-frame UV. The 
curve at z < 5 follows the He II history, which is gov¬ 
erned by quasars, with He II reionization completed at 
z ~ ‘i. The right panel compares /sfr, the critical SFR 
density from Equation (7), to the observed (and extrap¬ 
olated) SFR density. The increased LyC production ef¬ 
ficiency shifts the evolution in ionized volume fraction, 
Qi{z), toward higher redshifts. One can therefore offset 
QhyC ~ 0.006 against the combination of physical pa¬ 
rameters, (C'_f///esc)T' 4 ~° ®’^^ ~ 8.35 in this case. With 
the elevated LyC efficiency (model A), full reionization 
[Qi Ki 1) occurs at z « 7, and the “observed and ex¬ 
trapolated SFR exceeds /sfr at z < 9. The ionization 
history can be altered with various combinations of phys¬ 
ical parameters and by different extrapolations to low- 
luminosity galaxies. 

4.4. He I and He II Continuum Photon Luminosities 

From the same data used to determine LyC photon 
fluxes, we only need to change integration limits to com¬ 
pute fluxes in the He I and He II continua. The cal¬ 
culation of (5i(m), the He I-ionizing photon production 
for each initial mass, follows the same procedure as for 
Qo{m) using Equation (1) with Aum = 504.259 A. The 
Qi{m) curves with solar metallicity (Eigure 10) have a 
similar shape to the H I Lyman continuum curves in 
Figures 7 and 8. In both cases, the photon production 
increases with mass but flattens out at 120 Mq. Be¬ 
cause the photons are in the Wien limit {hiy ^ kT^s) of 
the Planck function, the intensity of radiation decreases 
exponentially in this short-wavelength range. Conse¬ 
quently the photon production rates for He I ionization 
are significantly less than the H I continuum rates, by up 
to an order of magnitude. The comparison between the 
new evolutionary tracks at two different metallicities fol¬ 
lows a different trend. The first difference is the behavior 
of the two most massive models in the sub-solar rotating 
model. Instead of flattening with increasing mass, the 
models seem to produce a rapidly increasing number of 
ionizing photons. The other sets of evolutionary tracks 
produce rather linear Qi{m) curves. 

From the data in Table 2, we find a factor of three range 
in Qnei depending on choice of IMF, metallicity, and ro¬ 
tation. For the “old” non-rotating models (Schaller et al. 
1992 tracks) we find QHei ~ (0.5 — 0.9) x 10®° photons 
per Mq at solar metallicity. At sub-solar metallicities, 
the He I continuum production is higher by 60% with 
QHei(subsolar)/Qnei(solar) = 1.6. In the newer non¬ 
rotating tracks (Ekstrom et al. 2012 and Georgy et al. 
2013) Qnei is a factor of 2 larger. For the new ro¬ 
tating models, the range is QHei ss (1.1 — 1.7) x 10®° 
photons per Mq at solar metallicity and 32% higher 
at sub-solar metallcity. Rotation elevates QHei at fixed 
metallicity for all three IMFs. At solar metallicity 
{Z = 0.014) the Ekstrom et al. (2012) tracks and corre¬ 
sponding atmospheres yield 94% higher efficiency, with 
(5Hei(i'ot)/(3Hei(nonrot) = 1.94. At sub-solar metallicity 
{Z = 0.002) the Georgy et al. (2013) tracks and cor¬ 
responding atmospheres were 22% more efficient, with 
(5Hei(i'ot)/QHei(nonrot) = 1.22. The overall He I effi- 




Fig. 10.— Lifetime-integrated He I continuum photon production 

in 10^^ photons vs. initial mass for several evolutionary 
tracks. Top: former standard tracks compared to new rotating 
and non-rotating solar metallicity tracks. “1992 Nonrot” tracks 
from Schaller et al. (1992), and “2012 Rot/Nonrot” tracks from 
Ekstrom et al. (2012). Bottom: new rotating and non-rotating 
tracks at solar and sub-solar metallicity. “2013 Rot/Nonrot” tracks 
are from Georgy et al. (2013). 


ciency ranges are comparable to those for QuyC- The 
recommended median value is Qnei = 0.6 x 10®° photons 
per Mq (non-rotating, solar metallicity) with significant 
enhancements (60-110%) at lower metallicity, unlike for 
the total LyC. The enhancements in Qiiei due to rota¬ 
tion are comparable to those for the LyC, at both solar 
metallicity (90% boosts) and sub-solar (22-35% boosts). 

Finally, we compute the He Il-ionizing photon fluxes 
Q 2 (jn) using Equation (1) with Aiim = 227.838 A. 
At higher ionization energies, the production of He II- 
continuum photons is several orders of magnitude lower 
than for H I or He I continua. Curves of (52(w) vs. 
initial mass are shown in Eigure 11 for rotating and non¬ 
rotating evolutionary tracks at solar and sub-solar metal¬ 
licity. Several features distinguish (52(w) from Qo{m) or 
Qi{m) curves. Eirst, the lower metallicity models pro¬ 
duce significantly more photons than models with higher 
metallicity. This is consistent with our expectations and 
the results of the Qo{m) and Qi{m) calculations. The 
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Fig. 11.— Lifetime-integrated He II continuum photon produc¬ 
tion (m) vs. initial mass. Top panel: Comparison of old 

evolutionary tracks (Schaller et al. 1992) and new solar metallic- 
ity rotating and non-rotating evolutionary tracks (Ekstrom et al. 
2012). Bottom: Comparison of new rotating and non-rotating evo¬ 
lutionary tracks at solar and sub-solar metallicity (Ekstrom et al. 
2012; Georgy et al. 2013). 


next distinction is that at sub-solar metallicity, the non¬ 
rotating models produce more He II continuum pho¬ 
tons than corresponding rotating models. This is dif¬ 
ferent from our expectations and the trends from solar- 
metallicity models. We can understand this discrepancy 
by examining the behavior of the sub-solar metallicity 
tracks, in which non-rotating models have a higher ef¬ 
fective temperature and luminosity for roughly 1.5 Myr, 
which constitutes a majority of their lifetimes. At the 
high energies of the He II continuum, the spectrum is 
quite sensitive to Teg. 

The recommended value for He II is Queii = I x 10®® 
photons per Mq (non-rotating, solar metallicity). The 
enhancements in QHeii due to rotation are a factor of 
1.48 at solar metallicity, but reduced by a factor of 
0.76 at sub-solar metallicity. Using the same pairwise 
comparisons as used above for QhyC and Quei, we find 
large boosts in Qiieii at lower metallicity. The ratio 
QHeii(subsolar)/(5Heii(solar) = 4.65 for old non-rotating 
tracks (and corresponding atmosphere models) as well 
as for the new rotating tracks. 


4.5. LyC Spectra in Model OB Associations 

The spectral evolution of an ensemble of stars of var¬ 
ious masses can be combined through an IMF and fol¬ 
lowed over time. Such spectral synthesis codes are valu¬ 
able tools for ultraviolet spectral libraries (Robert et al. 
1993; Rix et al. 2004; Leitherer et al. 2014) and stellar 
population synthesis codes such as Starburst99 (Lei¬ 
therer et al. 1999). Figure 12 shows a series of spectra for 
a cluster starburst containing 10® Mq in which the stars 
follow a Salpeter IMF (O.I < m < 120). The two panels 
show model spectra at times t = 0, I, 3, 5, and 7 Myr 
after a coeval burst of star formation, using evolution¬ 
ary tracks with rotation, at solar metallicity {Z = 0.014 
from Ekstrom et al. 2012) and at sub-solar metallicity 
{Z = 0.014 from Georgy et al. 2013). These models were 
created by Monte-Carlo sampling of stars for each mass 
for which we have an evolutionary track (m = 20, 25, 
32, 40, 60, 85, 120). For the Salpeter differential mass 
distribution, ^(m) = Km~°‘ with a = 2.35, the fraction 
of stars above mass m is given by 


N{>m) _ [m - rnmlx 


N 


[m 


-(a-l) 

min 


— (a—1)-| 

^max I 


( 8 ) 


The constant K is normalized to the total cluster mass 
M = 5.862Ar, and the total number of stars N = 
I6.582Ar for mass limits Wmin = O.I and rumax = _I20. 
For M = 10® Mq, we expect a mean number N = 
282,500 stars and mean stellar mass (m) = M/N = 
0.354 Mq. In many of our Monte-Carlo samples, the 
numbers fluctuate about this value, with small numbers 
at the high-mass end of the IMF. The mean numbers of 
stars at the high-mass end are: iV(> 60 Mq) r:! 30 and 
N{> 100 Mq) « 5. If the upper mass was extended to 
Wmax = 200, we would expect 40 stars above 60 Mq and 
15 stars above 100 Mq. As discussed earlier, the exis¬ 
tence of these very massive stars is controversial, owing 
to resolution effects. Moreover, the LyC from stars at 
m > 100 Mq may not escape the embedded cloud from 
which they formed, or from the dense gas produced in 
mass-loss episodes or binary mergers (Smith 2014). We 
therefore do not consider stars above the 120 Mq track. 

The spectra of each model are calculated from the 
star’s basic parameters on the evolutionary tracks, mul¬ 
tiplied by the number of stars for each mass and added 
to the total spectra. The stars are randomly chosen, 
weighted by the IMF, until the total mass of the clus¬ 
ter is 10® Mq. Owing to the rarity of 0-stars at the 
upper end of the IMF, some models had no stars above 
100 Mq. The LyC fluxes do not fade until after 5-7 Myr 
(Figure 13) owing to their post main-sequence evolution 
to O supergiants. However, fluxes at the shortest wave¬ 
lengths (A < 504 A) decrease after 5 Myr. 

As noted in Section 3.2, extragalactic observers have 
calibrated escaping LyC radiation from the flux-density 
ratios, F(9I2“)/F(9I2+) at the Lyman edge, and 
F(I500)/E(9I0), the far-UV-to-LyC ratio at 1500 A and 
900 A. These computed flux ratios are listed in Ta¬ 
ble 3 for model clusters calculated with rotating-star 
models at solar metallicity {Z = 0.014) and sub-solar 
metallicity {Z = 0.002). For coeval starburst mod¬ 
els, the intrinsic flux decrement at the 912 A Lyman 
edge ranges from F(9I2“)/F(9I2+) 0.4 — 0.7 over 
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Fig. 12.— Overlaid spectral distributions in specific luminos¬ 
ity, Lx in 10^® erg s“^ for Monte-Carlo models of a star¬ 

forming cluster. We assume 10^ Mq of stars with Salpeter IMF 
(0.1 < m < 120) and use evolutionary tracks with rotation and 
model atmospheres at the appropriate metallicities. Top: Solar 
metallicity {Z = 0.014, Ekstrom et al. 2012). Bottom: Sub-solar 
metallicity (Z = 0.002, Georgy et al. 2013). Each curve is offset by 
1.0 (in 10^® units) plotted, from bottom to top, at times £ = 0, 1, 
3, 5, and 7 Myr after a coeval burst of star formation. See Section 
4.5 for details. 


the first 5 Myr for both solar and sub-solar metallicities. 
The far-UV (1500 A) to LyC (900 A) ratio ranges from 
F(1500)/i^(900) 0.4 — 0.7 over the same time period. 

In their studies of escaping LyC from starburst galaxies 
at z ^ 3, Steidel et al. (2001) defined a relative fraction 
of escaping LyC photons normalized to escaping 1500 A 
photons: 

„ _ [T-(1500)/F(900)]int . . ... 

/esc,rel [i.(1500)/F(900)]obs ^ > 

where [r"(1500)/F(900)]int is the intrinsic ratio of non¬ 
ionizing, far-UV (1500 A) photons to ionizing photons 
(at 900 A) and riGM ,900 is the line-of-sight optical depth 
of the intergalactic medium to LyC photons at 900 A. 
The intrinsic spectral shape for massive stars above and 
below the Lyman edge is not well constrained by ob¬ 
servations, and one must rely on stellar atmosphere 
models and population synthesis. Both Steidel et al. 



Fig. 13.— Monte-Carlo simulations, as in Fig. 12, of the time 
evolution of LyC production rate Qo{t) [photons s“^] in a clus¬ 
ter with 10^ Mq, formed coevally in a Salpeter mass distribution 
(0.1 < m < 120). We use evolutionary tracks at solar metal¬ 
licity with rotation (Z = 0.014, Ekstrom et al. 2012) and show 
two runs with different selections of stellar masses, from initial 
random-number seeds. Late-time spikes (3.5 — 6.2 Myr) arise from 
the post-main-sequence evolution of stars with 40-60 Mq. 


(2001) and Inoue et al. (2005) adopted an intrinsic ratio 
F(1500)/U(910) = 3.0. Because our (1500 A/900 A) ra¬ 
tios are considerably lower (0.4-0.5) this will reduce the 
inferred escape fractions in high-redshift starburst galax¬ 
ies, which are proportional to the assumed intrinsic ra¬ 
tio. This decrease in the intrinsic ratio may resolve the 
paradoxical inference (Shapley et al. 2006), who found 
a corrected flux ratio, U(1500)/F(900) ~ 2.9 implying 
that /esc is greater than 100%. After 5-7 Myr, when the 
most massive 0-type stars evolve off the main sequence 
and eventually die, the (1500-to-900 A) flux ratios begin 
to increase. 

5. DISCUSSION AND CONCLUSIONS 

This work provides an update of the production rates 
of ionizing photons from massive stars, based on a new 
generation of stellar evolutionary tracks coupled to mod¬ 
ern stellar atmosphere models. For a range of parameters 
(tracks, metallicities, rotation), we find a 50% increase in 
the LyC production efficiencies compared to previous cal¬ 
ibrations. The new stellar tracks from the Geneva group 
were calculated for both solar and sub-solar metallici¬ 
ties, and they also examined the effects of stellar rota¬ 
tion with surface velocities 40% of breakup. For the grid 
points (reff,log 5 ) on these tracks, we computed model 
atmospheres with the non-LTE code WM-basic, which 
treats 3D expanding atmospheres and line-blanketing by 
lines of heavy elements whose metallicities were consis¬ 
tent with those in the stellar tracks. We computed the 
production rates of ionizing photons in the ionizing con¬ 
tinuum of hydrogen, as well as of He I and He II, and their 
dependence on metallicity, rotation, and initial stellar 
mass. Next, we integrated these rates over three different 
IMFs (Salpeter, Kroupa, Chabrier) to derive LyC photon 
production efficiencies, QuyCi in units of LyC photons 
per solar mass of star formation. 

Our calculations used the recent grids of evolution¬ 
ary tracks from the Geneva group, who made specific 
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assumptions about stellar mass loss rates and rotation. 
These grids could be revised in the future, to account 
for improved parameterization of mass-loss rates, binary 
evolution, and the upper extent of the IMF. Recent ob¬ 
servations of massive stars have introduced uncertain¬ 
ties in mass-loss rates, particularly in late (WR) stages 
(Grafener & Vink 2013; Smith 2014). Binaries could af¬ 
fect the stellar populations owing to tides, mass transfer, 
and mergers (de Mink et al. 2013). Some effects of stellar 
mergers are effectively included in tracks with increased 
rotation. After considering a range of metallicities and 
the effects of stellar rotation, the main conclusions of our 
study can be summarized as follows: 

1. By coupling new stellar evolutionary tracks at sev¬ 
eral metallicities, with and without rotation, to 
non-LTE model atmospheres, we find a range of 
LyC photon production efficiencies QhyC = [3.1 — 
9.4] X 10®° photons per Mq of star formation. The 
median LyC production efficiency is (6 ± 2) x 10®° 
LyC photons per Mq, an average increase of 50% 
from previous models. 

2. The new LyC efficiency is equivalent to a 
SFR calibration of 10®°°=*=°'2 photons s“^ per 
Mq yr“^. With these new rates, the criti¬ 
cal SFR density in a clumpy medium can be 
written psFR = (0.012 Mq yr“^ Mpc“®)[(l -|- 

z)/8]3(CH/3)(0.2//e,e)(0.006/QLyc)T4-°•®4^ 
scaled to fiducial values of IGM clumping factor 
Ch = 3, LyC escape fraction /esc = 0.2, and 
electron temperature Tg = lO"* K. 

3. The higher LyC efficiencies make it easier to com¬ 
plete reionization by z w 7 for various physi¬ 
cal parameters Ch, /esc, 7/ and extrapolations 
of the galaxy luminosity function to magnitudes 
well below those observed in the HST/WFC3 deep 
fields. For LyC efficiency QuyC = 6 x 10®° pho¬ 
tons per Mq, reionization is complete by z = 7 for 
{Ch/Usc) « 15 - 30 and = (1 - 2) x 10^ K. 

4. We also computed efficiencies for He I and He II 
continua. For He I, we find Quei = 0.6 x 10®° 
photons per Mq (non-rotating, solar metallicity) 
with substantial (60-110%) enhancements at sub¬ 
solar metallicity and with rotation (90% boosts at 
solar metallicity and 22-35% boosts at sub-solar). 
For He II, we find Qneii = 1 x 10®® photons per 
Mq (non-rotating, solar metallicity) and a factor 
of 4.65 higher at low metallicity. 

5. Because of the H I opacity in the upper stellar at¬ 
mosphere, non-LTE models of stars with < 
50, 000 K exhibit 40-60% intrinsic LyC edges at 
912 A. These should be visible with next-generation 
UV missions that can observe escaping radiation in 
the rest-frame LyC of low-z starburst galaxies. 

6 . Monte-Carlo simulations of starbursts find ratios 
of far-UV (1500 A) to LyC (900 A) flux density 
E(1500)/F(900) w 0.4—0.5 over the first 5 Myr, for 
models with metallicity Z = 0.014 and 0.002. The 
flux ratio at the Lyman edge E(912“)/E(912+) w 
0.4 — 0.7 over the same time interval. 


Several improvements could be made in future calcula¬ 
tions. Additional stellar tracks at initial masses between 
25 and 60 Mq, with rotation and super-solar metallicity 
would strengthen our results for the metallicity depen¬ 
dence of QhyC- They would also provide a further basis 
of comparison to older super-solar evolutionary tracks 
(Schaller et al. 1992). Models with a range of rotation 
rates would be helpful, since rotation speeds at 40% of 
breakup may be excessive and produce ionizing SEDs 
inconsistent with nebular emission lines (Levesque et al. 
2012). Our calculations rely on rotating tracks in which 
all stars have the same initial rotation. None of the 
current tracks include a description of magnetic fields, 
although that has not yet been shown to have a signifi¬ 
cant effect. Future stellar population models will treat 
rotation more realistically, with better understanding of 
their ionizing photon budgets. 

Several changes in our computational scheme could 
also result in more accurate photon fluxes and pro¬ 
duction rates. When creating the grid (log 5 , Teff) of 
model atmospheres, we continued adding models until 
we reached convergence in photon production rates. 
Adding more model atmospheres, up to one per time 
step of the evolutionary tracks, would create more preci¬ 
sion. Another small effect is the contribution of ionizing 
photon production of lower mass stars. We currently 
consider models with masses m > 25 Mq. Lower mass 
stars produce an order of magnitude or more fewer LyC 
photons. Filling in the grid at other metallicities and 
for additional masses between 35 M 0 and 120 Mq would 
be helpful. For reionization, the most important factors 
that could increase the total LyC production are the 
stellar populations at low metallicity, the extent of the 
upper main sequence, the influence of rapid rotation, 
and the role of binaries. 
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TABLE 1 

Ionizing Photon Production^ vs. Initial Stellar Mass 


Model Tracks*’ 

(Z/Zq) 

Qo(120 Mq) 
(10®3 phot) 

Qo(85 Mq) 
(10®3 phot) 

Qo(60 Mq) 
(10®3 phot) 

Qo(40 Mq) 
(10®3 phot) 

Qo(32 Mq) 
(10®3 phot) 

Qo(25 Mq) 
(10®3 phot) 

Ekstrom (Rot) 

0.014 

11.5 

9.34 

6.32 

3.19 

2.06 

1.03 

Ekstrom (Non-Rot) 

0.014 

7.49 

5.11 

3.58 

1.58 

0.902 

0.445 

Georgy (Rot) 

0.004 

16.4 

9.95 

5.70 

2.59 

1.76 

1.67 

Georgy (Non-Rot) 

0.004 

9.95 

7.31 

4.47 

2.19 

1.43 

0.871 

Schaller (Non-Rot) 

0.020 

7.91 

5.78 

3.66 

2.42 


0.474 

Schaller (Non-Rot) 

0.004 

8.29 

6.46 

4.47 

2.18 


0.648 


^ Total number of LyC photons, Qo(^) iR 10®^ photons, produced over stellar lifetime for stars of initial mass m. Evolutionary tracks and 
non-LTE model atmospheres are described in Section 3 and Table 2 footnote. 

^ Evolutionary track references: Schaller et al. (1992); Ekstrom et al. (2012); Georgy et al. (2013), of various metallicities (Z/Zq) and 
rotating or non-rotating models as noted. 


TABLE 2 

Ionizing Photons^ with Various IMFs, Tracks, and Metallicities (Z) 


Model 

IMF*’ 

Tracks*’ 

Rot-Type 

'Rimiii 

(Me) 

R^max 

(Me) 

Slope 

a 

Metals 

Z 

QlvC 

(10^1) 

^Hel 

(10^0) 

QHell 

(10®®) 

Solar 

1 

4 

Nonrot 

0.1 


2.35 


0.367 

0.552 

0.918 

Solar 

1 

5 

Rot 

0.1 


2.35 

0.014 

0.603 

1.11 

1.73 

Solar 

1 

5 

Nonrot 

0.1 

120 

2.35 


0.322 

0.571 

1.17 

Solar 

2 

4 

Nonrot 


120 

2.3 

0.020 

0372 

0.864 


Solar 

2 

5 

Rot 


120 

2.3 

0.014 

0.940 

1.74 


Solar 

2 

5 

Nonrot 


120 

2.3 

0.014 

0.503 

0.895 


Solar 

3 

4 

Nonrot 



2.3 


0.370 

0.558 

0.935 

Solar 

3 

5 

Rot 



2.3 

0.014 

0.607 

1.13 

1.78 

Solar 

3 

5 

Nonrot 



2.3 


0.325 

0.578 

1.21 

Sub-solar 

1 

4 

Nonrot 

0.1 

120 

2.35 

0.004 

0.399 

0.884 

4.25 

Sub-solar 

1 

6 

Rot 

0.1 

120 

2.35 

0.002 

0.583 

1.47 


Sub-solar 

1 

6 

Nonrot 

0.1 

120 

2.35 

0.002 

0.433 

1.21 


Sub-solar 

2 

4 

Nonrot 


120 

2.3 

0.004 

0.623 

1.38 

6.75 

Sub-solar 

2 

6 

Rot 


120 

2.3 

0.002 


2.31 

12.8 

Sub-solar 

2 

6 

Nonrot 


120 

2.3 

0.002 



16.9 

Sub-solar 

3 

4 

Nonrot 


120 

2.3 




4.36 

Sub-solar 

3 

6 

Rot 


120 

2.3 



1.49 

8.28 

Sub-solar 

3 

6 

Nonrot 


120 

2.3 


0.436 

1.22 



^ Photon production,QLyC, Qneii QHeii iR iRst three columns are for ionizing continua of H I, He I, and He H, respectively. Total LyC 
production (Qlyc) is expressed in units of 10®^ lifetime-integrated photons per Mq of star formation, calculated for various IMFs and 
metallicities (columns 1-8) and for evolutionary tracks with or without rotation. Values of QhcI RRd QhcII RI"® iR RRits of 10^^ and 10^^ 
photons per Mq, as noted in header. 

^ IMF references: (1) Salpeter (1955); (2) Kroupa (2001); (3) Chabrier (2003). Evolutionary track references: (4) Schaller et al. (1992); (5) 
Ekstrom et al. (2012); (6) Georgy et al. (2013). Using the Kroupa or Chabrier IMFs increases the LyC production relative to the Salpeter 
IMF by factors of 1.56 and 1.64, respectively. 


TABLE 3 

Flux Ratios in Composite Cluster Spectra^ 


Time® 

F(912-)/F(912+) 

Z = 0.014 

F(912-)/F(912+) 

Z = 0.002 

F(1500)/F(900) 

Z = 0.014 

F(1500)/F(900) 

Z = 0.002 

0 Myr 

0.54 

0.58 

0.472 

0.455 

1 Myr 

0.57 

0.71 

0.457 

0.438 

3 Myr 

0.65 

0.71 

0.506 

0.404 

5 Myr 

0.39 

0.41 

0.777 

0.727 

7 Myr 

0.50 

0.31 

0.630 

1.03 


^ Ratios of flux densities, (erg cm ^ s ^ A ^) at key wavelengths: E(912 )/E(912'*") is the flux decrement at the 912 A Lyman edge; 
E(1500)/E(900) is the ratio of fluxes in the far UV (1500 A) and 900 A, just shortward of the LyC edge (Steidel et al. 2001; Shapley et al. 
2006). These ratios are shown for simulated clusters (Figure 12) with solar and sub-solar metallicities, as noted, using evolutionary tracks 
with rotation (Ekstrom et al. 2012; Georgy et al. 2013). 

^ Time in Myr following a coeval burst of star formation in a cluster with 10^ Mq with a Salpeter IMF (0.1 < m < 120) as shown in 
Figure 12. 






































